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ABSTRACT 

We describe a simple method for simulating the dynamics of small grains in a dusty gas, 
relevant to micron-sized grains in the interstellar medium and grains of centimetre size and 
smaller in protoplanetary discs. The method involves solving one extra diffusion equation 
for the dust fraction in addition to the usual equations of hydrodynamics. This “diffusion 
approximation for dust” is valid when the dust stopping time is smaller than the computational 
timestep. We present a numerical implementation using Smoothed Particle Hydrodynamics 
(SPH) that is conservative, accurate and fast. It does not require any implicit timestepping and 
can be straightforwardly ported into existing 3D codes. 

Key words: hydrodynamics — methods: numerical — protoplanetary discs — (ISM:) dust, 
extinction — ISM: kinematics and dynamics 


1 INTRODUCTION 

Small grains rule the interstellar medium (ISM). Micron-sized dust 
grains absorb ultraviolet radiation from hot, young stars and re-emit 
it in the infrared. Understanding how these grains interact with the 
gas is critical to understanding both the dynamics and thermody¬ 
namics of the ISM, and to interpreting observational results which 
usually assume a fixed gas-to-dust ratio in order to derive physical 
quantities such as the gas column density. 

Modelling such grains presents a severe computational chal¬ 
lenge, since small grains are tightly coupled to the gas by the mu¬ 
tual drag force. This presents both a short timescale problem, since 
the stopping time of the grains is much shorter than the typical com¬ 
putational time, and a short lengthscale problem, since the physical 
separation between the dust grain population and the gas is much 
smaller than typical distances in the ISM. _ 

In a recent series of papers (iLaibe & Pricel2012'allRl2014all5l3) 
we have outlined the limitations associated with modelling dust and 
gas using the standard two fluid approach, where they are regarded 
as separate fluids coupled by a drag term. Typically the gas is rep¬ 
resented by a set of particles or grid cells, while the dust is rep¬ 
resented by a separate set of pressure-less particles coupled to the 
gas by a drag term. The length and timescale problems discussed 
above mean that with this approach one needs both infinite spa¬ 
tial and temporal resolution to accurately ca pture the dynamics of 
small g rains in the limit of perfect co upling Juaibe & Pricell201^ : 
but see iLoren- Aguilar & Batel (l20l4) for an alternative approach). 
However, this is the limit in which the mixture can be accurately 
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described as a single fluid movin g at the baryce ntric velocity. In 
iLaibe & Pried (1201 4alf9) f hereafter lLP 1 4al : IlP 1 4bh we showed how 
the equations for a coupled dust-gas system can be reformulated to 
describe this single fluid mixture without loss of generality, solv¬ 
ing both the length and timescale issues and also preventing ar- 
tiflcial trapping of dust pa rticles below the resolution of the gas 
(lLP12al : lAvliffe et al.ll2012h . The method is similar to the approach 
to other multi-fluid syste ms in astrophysics such as ionised plasmas 
JPandev & Wardlell2008h . but more general since it can be imple¬ 
mented without any approximations. 

In lLPlil we derived a Smoothed Particle Hydrodynamics 
(SPH) algorithm based on the fully general one fluid method and 
showed that it could accurately capture the dynamics of dust-gas 
mixtures in both the weakly coupled and tightly coupled limits. For 
problems involving small grains, however, the full machinery of the 
one fluid formulation is unnecessary and a much simpler and com- 
putationall v inexp ensive approach is possible, as outlined in Sec¬ 
tion 3.3 of lLP14aL This approximation is accurate when the stop - 
ping time, 4, is less than the Courant timestep (Eq. 115 in lLP14ah . 

Our goal in this paper is to derive a numerical implementation 
of this much simpler formulation, since there are many situations 
in astrophysics where the dynamics of small grains is the dominant 
effect. This includes simulations of galaxies, star formation in the 
interstellar medium — where small grains control the thermody¬ 
namics — and the settling and migration of dust in protoplanetary 
discs. We summarise the analytic formulation and its applicability 
in Sec. 121 the numerical implementation is described in Sec. [3] and 
tests are presented in Sec.|4] A public version of the ndspmhd code 
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2 Price & Laibe 


(v2.1) implementing the algorithms and with the precise setup of 
the test problems is released alongside this papefl 


2 THE DIFFUSION APPROXIMATION FOR DUST 


2.1 Continuum equations 


2.1.1 General case 


In iLPUal we showed that, to first order in 4/r, where T is the 
timescale for a sound wave to propagate over a typical distance L, 
the equations describing the evolution of a dust-gas mixture can be 
written in the form 


dt 

dv 

dt 

&u 

d? 


— = -p(V. V), 


t: - (1 “ ^)/g + ^/d + /» 


--V • [6(1 - e)pt,^f ], 

P 

P 

(V • v) -|- 6t^ (Ay* • V) U + Aheat “ Acoob 
Pg 


( 1 ) 

( 2 ) 

(3) 

(4) 


where p is the total density of the mixture, e = p^/p is the mass 
fraction of dust, f represents accelerations acting on both compo¬ 
nents of the fiuid while /g and /d represent the accelerations acting 
on the gas and dust components, respectively. A/ = /d - /g is the 
differential acceleration between the gas and dust, u is the specific 
thermal energy of the gas, P is the gas pressure, and Aheat and Acooi 
are additional heating and cooling terms, respectivel}|l . The veloc¬ 
ity V is the barycentric velocity of the mixture, defined as 


^ ^ PdVd+PgVg 


: 6Vd + (1 - 6)Vg 


(5) 


In the so-calle d _ ter minal velocit y approx i mation 

(lYoudin & GoodmanI 20051: IChiand l2008l : |b arranc 3 l2009l: 
iLee et al.l20ld : ljacauet et al.11201 ih assumed in Equations [T}01 A/ 
is rapidly balanced by the drag. Thus, the time dependence of the 
differential velocity can be ignored, and the differential velocity 
between the gas and dust is given by 


Av = (Vd - Vg) ^ tsAf. (6) 

This also implies that the anisotropic pressure term in the momen¬ 
tum equation (see lLP14^ should be neglected. The terminal veloc¬ 
ity approximation is valid when the drag coefficient K is large such 
that the stopping time, 

^ PdPg _ 6(1 - 6)p 
^ 7r(pd+Pg) K ’ 

is short compared to the timestep. Various physical prescrip- 
tions for K in the Ep stein and Stokes drag regimes are given in 
iLaibe & Pried (12012bh but the essential point is that K is inversely 
proportional to the grain size, being large for small grains. 

The differential acceleration A/ depends on the physics in the 
problem, i.e. the forces affecting the gas but not the dust, which may 


^ http://users.monash.edu.au/~dprice/ndspmhd/ 

^ Eq. |4| di ffers fro m the expression we gave for the “first order approxi¬ 
mation” in lLP14aL The drag heating term, eAv^ I ts, is clearly negligible in 
the terminal velocity approximation and the PdV work term should involve 
V • V rather than V • Vg. Both approximations are required for the numerical 
scheme to conserve total energy as defined in the terminal velocity approx¬ 
imation (Eq.[^. 


include pressure, magnetic and other forces. In our numerical im¬ 
plementation we consider the contributions from the pressure gra¬ 
dient (see below) and also the artificial viscosity term, which should 
likewise affect the gas only. 


2.1.2 Hydrodynamics 


For the simple case of hydrodynamics, the only force is 
gradient, giving 

the pressure 

VP 

/g = -—; /d = o, 

Pg 

and thus 

(8) 

VP 

A/=—, 

Pg 

(9) 

giving Equations [TJS] in the form 

1 = -P(V-v), 

(10) 

dv VP 

dt p 

(11) 

J = -Iv-(e4VP), 

dt p 

(12) 

du P ep 

= - (V • v) H- (VP • Vw) -1-Aheat “ Acooi- 

(13) 


dt pg pg 


These are similar to the usual equations of hydrodynamics in the 
absence of dust. The only differences are the extra equation that 
describes the evolution of the dust fraction; the modifications to the 
thermal energy equation; and the fact that the pressure is related 
to the gas density only, not the total density (see Sec. [2X3] below; 
this giv es the zeroth order effect of a ‘heavy fluid’, as discussed in 
iLPUah . 


2.1.3 Equation of state 

The equation set is closed by the usual equation of state specify¬ 
ing the gas pressure P in terms of the gas density and temperature. 
Unless otherwise specified in this paper we assume an adiabatic 
equation of state, i.e. 

P = {y- l)pgu = (y- 1)(1 - 6)pw, (14) 

where y is the usual adiabatic constant. 


2.2 Timestepping 

The main change when adopting the formulation given above com¬ 
pared to hydrodynamics is the addition of the diffusion equation for 
the dust fraction 03 . This introduces an additional constraint on 
the timestep when the diffusion coefficient is large. Assuming an 
isothermal equation of state P = c^pg = cl(l - e)p and a constant 
density, can be written as a simple diffusion equation for e 

^=V-(rjVe), (15) 

dt 

where the diffusion coefficient q = et^c]. This implies a stability 
constraint of the form 

hf 

At < At^ = Co — = Co- T, (16) 

q etscj 
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A fast algorithm for small grains in SPH 3 


where Cq is a dimensionless safety factor of order unity and h is the 
resolution length (the smoothing length in SPH). We can rewrite 
([T^ as 

(17) 

where C is a constant and Atcour = Coh/cs is the usual Courant 
condition. This implies that the timestep is constrained when the 
stopping time is long — the opposite of the usual situation where 
the timestep is constrained when the stopping time is short. This is 
the main advantage of using the diffusion approximation — small 
grains can be integrated explicitly. 

Specifically, the diffusion timestep becomes the limiting 
timestep when 

6ts > Atcour- (18) 

However, this is also the criterion for whe n the terminal veloc¬ 
ity approximation breaks down fsee lLPljj) . This implies that the 
diffusion approximation becomes inaccurate precisely when the 
timestep implied by ([Tfil) starts to constrain the timestep, because 
at this point the time-dependence in Av becomes important. Once 
thi s occur s, one should revert to the general formulation given 
by lLP14bl where Av is explicitly evolved, or a two-fluid method. 
Physically this transition occurs once grains grow beyond a certain 
size, implying that the stopping time becomes long, or equivalently 
when one has enough temporal resolution to resolve the timescale 
on which the differential velocity is changing. 


2.3 Validity of the diffusion approximation for astrophysics 


Under what circumstances is the diffusion approximation valid for 
astrophysics? Consider a drag force described by the linear Epstein 
regime, appropriate to small grains at l ow Mac h number. In this 
case the drag coefficient is given by (e.g. lLP12bl) 


K 


— PgPd 


47r “^grain 

3 Wlgj-ain 



(19) 


where ^grain is the grain size and mgrain is the mass of an individual 
grain. Assuming mgrain = f^Pgrain^^gram’ where pgrain is the intrinsic 
grain density, the stopping time is 


4 — 


Pgrain “pgrain 
pCs 



( 20 ) 


2.3.1 Grains in the interstellar medium 

Evaluating this for dust grains in a molecular cloud, we have 


4 = 2.5xl0^yr 


Pgrain 

Ig/cm"^ / \ 0.1pm / \ 10"^^g/cm' 


0.2km/s 


( 21 ) 

This indicates that the diffusion approximation is valid for small 
grains in the interstellar medium, since the stopping time is much 
smaller than the dynamical time (~ lO^yr). 


2.3.2 Protoplanetary discs 


Eor a protoplanetary disc, the relevant comparison is to the orbital 
timescale since the pressure timescale Hjc^ = 1/Q. A reasonable 
criterion for validity is therefore that 




Pgrain “pgrain 

E 


1 . 


( 22 ) 


This suggests the approximation is valid for grain sizes 


“pgrain 


lO^cm 


E 

lO^g/cm^ 


Pgrain 

Ig/cm^ 



(23) 


Hence diffusion is a reasonable approximation for grains of ~cm 
size and smaller in protoplanetary discs. This maximum size is 
smaller in the outer disc regions, since typically the surface den¬ 
sity is inversely proportional to distance from the central star. We 
examine this experimentally in Section 1441 


3 IMPLEMENTATION IN SMOOTHED PARTICLE 
HYDRODYNAMICS 


3.1 Implementation using two first derivatives 

The SPH repres entation of a more general form of Eqs.[T}|4|have 
been derived in |LP14bl and so our first approach is to adopt the 
same discretisation but with Av prescribed by Eq. i giving 


dVg 

dt 

dCg 

dt 


dUg 

dt 


^ ^b^ab(ha)-> 


(24) 


b 

+fa^ 

b 


1 


Pa+^fla 


^a(l “ ^a)is,a 


Pb + 

-^aWgbQlg) + --^VgWgtih) 


^aPa 


^hPl 

A/. • '^aWabiha) 

Afb ■ VaW^bihb) 


(25) 


^a is,a 

a 


ClbPb 

^ nibiPa + cfla) {Va ” Vfc) • V„W„fc(/l„) 


(26) 


^aPa 


a/«E f^bi,^a ^b^^ a^dlg}j{hg)^ 


(27) 


where Wg^ is the usual SPH kernel (we use the usual cubic spline 
kernel throughout this paper unless otherwise indicated), h is the 
smoothing length, Q is the usual term related to smoothing length 
gradients 


Og^l- 


dhg ^ dWgl,{hg) 

dpg ^ dhg 


(28) 


and h is related to p in the usual manner requiring an iterative pro¬ 
cedure to solve EQ. I^(lLP14bl:IPrice & Monaghanil2004l2007h and 
unl ess otherwis e specified we use a ratio of h to particle spacing of 
1.2 (|Pricell2012h . The reader will notice that the first two equations 
are identical to the usual density summation and momentum equa¬ 
tion in SPH. The only differences, mirroring the continuum case 
fEas. fTMT^ . are the addition of the diffusion equation m for the 
dust fraction, the extra terms in the thermal energy equation (IZTI ) 
and the dependence of the pressure on the gas density rather than 
the total density in the equation of state ([T4]) . 

The differential force between the fluids implied by our for¬ 
mulation of EQ.l25lis 


A/g = (29) 

where 


(i-f„)/“ = -E 


mt 




ab,a 


^aPl 


^aWabiha) + 




ab,b 


^bPi, 


VaWabihb) 

(30) 
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4 Price & Laibe 


This A/, computed as above, is then used to evaluate Equations [26l 
and [571 requiring a separate loop over the particles. 


3.2 Shock-capturing terms 

3.2.1 Artificial viscosity 

We formulate the artifici al visco sity term following the more gen¬ 
eral algorithm derived in lLP14bl but slightly modified to appear as 
separate qa and qt terms to avoid averaging the kernel gradients, 
following the formulation of artificial viscosity used in t he Phan¬ 
tom code (jPrice & Federrathir201Ql : lLodato & Pricell201Qh . We use 


^ab,a 


2Pa(l ^a)^sig,a^ab ' ^ab- 

0 


^ab • TaZ; < 0 
^ab ’ ^ab — 0 


(31) 


where Vab = Va - Vi, (similarly for Tab) and the signal speed Vsig 
corresponds to the usual choice for hydrodynamics, i.e. 

l^sig,a “ ^aCs,a (3\^ab ' ^ab\’ (32) 

where a e [0,1] is the linear dimensionless viscosity parameter 
(in general this can be individ ua l to each particle, e.g. w hen using 
the lMorris & Monaghan!! 19971 or ICullen & Dehnenll201Ql switches) 
and p (typically = 2) is the Von Neumann-Richtmeyer viscosity 
parameter. 

The q^^ term and the signal speed invol ve the ju mp in total 
velocity rather than the gas velocity, unlike in lLP14hl where only 
the gas velocity is used. This is both physical and practical: In the 
terminal velocity approximation the difference 


V - Vg = 6tsA/, (33) 

is small by definition. The practical side is that it we do not know 
A/ prior to the evaluation of Eq. [251 so it is not possible to use 
the gas velocity directly in the artificial viscosity term without an 
iterative approach. 


3.2.2 Artificial conductivity 

We write the artificial conductivity term, neces sary for correct treat- 
ment o f contact discontinuities dPricel [200^ . similar to that in 
lLP14hl giving 


(dua\ 

\ /cond 


1 


1 - fa 




Fabiha) 




Qab,b 

^bpl 


Fabifib) 


(34) 


where = Fab^ab and 

Qab,a — 2 ^MP«^sig,u(^a “ ^Z?)? (35) 


with au e [0,1 ] the dimensionless conductivit y parameter and 
Vsig,u = \^ab • ^ab\ (iPricell2008l : IWadslev et al.ll2008h . 


3.3 Conservation properties 

Equation [24] manifestly conserves the total mass since the mass of 
the SPH particles is constant. Similarly it can be straightforwardly 
verified that the total momentum is conserved, since 


due to the fact that the resulting double summation is antisymmetric 
in the particle indices a and b. Likewise the total angular momen¬ 
tum is conserved, since 


d v-n v-n dVa ^ 

^ / 1 ~ / ^ ^aXa ^ “ 0. 


(37) 


dt^ “ “ - 

a a 

(for more details, see Equation 33 in !PricelEoi2h . Finally, one may 
also verify that the total mass of each species is conserved, since 

a 

The proof is identical to that given in !LP14bl and again results from 
the fact that the double summation is antisymmetric with respect to 
the particle indices. 

The total energy o f the mi xture in the terminal velocity ap¬ 
proximation is given by ZLPUah 


J + Pgwj f + P(1 - e)M 


dV. 


(39) 


This is simpler than the full one fluid expression (Eq. 61 in lLP14a|l 
as the term involving Av^ can be neglected. Discretised onto the 
mixture particles, the energy becomes 




1 


-Vl + {l- €a)Ua 


Conservation of energy implies that 


dE 


dt 


= E' 


dVa ,, ,dMa dea 


= 0 . 


(40) 


(41) 


Substituting Equations!^ andin the above, we require for en¬ 
ergy conservation that 




dt 


-EE 

a b 

-EE 


nib 


niamb 


maMb 


P + a^^ 

^ ^ab,a ^ . 

2 ^fid/abO^a) 


FlaPl 


ab,b 


^hp\ 

^a(l “ ^a)^fis,a 


-EE niamb 

a b 


^aPa 

^a(l “ ^b^^bis,b 


Va • '^aWabifib) 

^fa • ^aWabiha) 


ClbPb 


^fb • ^oWabQlb) 


Swapping the summation indices a and b in the second and fourth 
terms, using the antisymmetry of the kernel gradient VbWbaiha) = 
-VaWabiha) and collecting terms we have 

V PI V V f^«"^^aZ^ 

2^ mail - ^ 2^ 2^ maMb 


■EE"*- 




ap^ 

(1 “ fa)fa4,( 


F^aPa 


(Va - Vb) • ^fid/abiha) 

{Uq ~ ^b^dsfa • ^fi^abFbfi) 


(42) 


from which it is straightforward to verify that, with duajdt given by 
Eq.!27] total energy is conserved exactly. 

Thus, the approximate version of the one fluid algorithm re¬ 
tains all of the conservation properties of both the original SPH 
method and the general one fiuid approach derived in !LP14bl 


3.4 Implementation using direct second derivatives 




The main disadvantage of the formulation given above is that it re¬ 
quires a third loop over the particles to compute the de/dt term. 
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beyond the two loops required for the density and force, respec¬ 
tively. This is because A/ is required before Eq. can be eval¬ 
uated, but must be computed after the right hand side of (|23 is 
known. Thus in general this scheme is 1/3 more expensive than 
a standard SPH code. Here we provide an alternative scheme that 
does not require this extra loop. The two implementations are com¬ 
pared in Section m 


3.4.1 Dijfusion equation for the dust fraction 

We can avoid the extra loop over the particles by discretising the 
second derivative in Eq. [3] directly, similar to the usual way that 
dissipative terms are treated in SPH. To do this we assume that 
viscous forces do not significantly drive the differential velocity 
between the fluids, i.e. that A/ is given by Eq.[9]and therefore that 
Eq. [3] is given by Eg. \\2\ We then discretise Eq. [12] in the usual 
manner following IClearv & Monagh^ (Il999h : 


dt 



(43) 


where D = 6^, Fgb = \[Fab{ha) + Fabih)] and Fgb is defined 
such that '^Wab = Fgb^ab- It is straightforward to show that this 
expression also conserves both the total mass of dust and gas, since 
the resulting double summation in Eq. [38] is antisymmetric with 
respect to the particle index. 


3.4.2 Harmonic V5. arithmetic mean 

In the original IClearv & Monagha^ (Il999h paper (see also 
iMonaghanI l2005h it was suggested to use the harmonic mean in¬ 
stead of the arithmetic mean of the diffusion coefficient, i.e. 


3.4.3 Thermal energy equation 


In order to conserve energy, the corresponding expression for dujdt 
when using Equation |43] for dejdt is given by 


dUg 

dt 


——I- — y mbiPa + cflf (»'a - »'fc) • '^aWabiha) 

^a)Pg ^ 

^^ — (Wa “ Ub){F>g Df){Pg - P^)—^45) 

2(1 - eg)Pg ^ Pb \rgb\ 


At first sight the second term is a rather strange one and it is not 
at all clear that this should translate to the correct physical term in 
Eq.|27] Yet, amazingly, it does — the proof is given in Appendix lAl 
Hence there is no disadvantage in using this alternative formulation 
with respect to conservation properties. The shock capturing terms 
remain the same as in Section [T2] 


3.4.4 Choice of smoothing kernel 

Although the formulation of second derivatives in SP H using the 
kerne l gradient P3]) is now more than 30 years old (iBrookshawl 
Il985h . and while it is clearly better than using directly, to our 
knowledge there has been no systematic investigation of the best 
kernel to use in order to compute a second derivative. In particu¬ 
lar, on the dust settling test in Section |4!4] we found that using (l43l ) 
with the cubic spline could give quite noisy results. Hence for this 
test we instead adopted the quintic kernel instead (see Sec. HU). 
While this results in a more accurate estimate, it is also more ex¬ 
pensive due to the larger kernel radius. Hence a more systematic in¬ 
vestigation of suitable kernels for s econd d erivatives in SPH would 
be valuable here. Eor example, in lLP12al we found double-hump 
shaped kernels to be an order of magnitude more accurate com¬ 
pared to standard kernels for computing the drag terms in the two 
fluid method at no additional cost. 


3.4.5 Two first derivatives v>s. direct second derivatives 


^ _ V (p _ p \ 

dt ^ pgPb {Dg + Db) ^ ^ \rgb\ 


(44) 


with the motivation being that this better handles the case where the 
diffusion coefficient D is discontinuous. However, we found this 
could give incorrect results. Imagine the dust confined to a layer 
such that €g = 0 for some particle, a, outside the layer, with 6^ ^ 0 
for particles inside the layer. In this case the harmonic mean is zero 
for every pair involving particle a since dcg/dt is always zero. Thus 
it is impossible for the layer to move into the region where e was 
initially zero, which is clearly incorrect (consider for example a dis¬ 
crete layer of dust descending under gravity). With the arithmetic 
mean we find no such problem and it is easy to prove that the for¬ 
mulation is correcl^ for example with a procedure similar to the 
one we use in Appendix O 


To our knowledge there exists no systematic study on whether 
it is better to compute second derivatives in SPH directly or us¬ 
ing t wo consecutive first derivatives (though see I Watkins et al.l 
1 19961) . In principle both approaches yield a second order ap¬ 
proximation provided that the particles are well ordered, and in 
the context of implementing physical visco sity terms in SPH 
both approaches ha v e been advocated ( e .g. lElebbeetak 1 1994 ; 
Watldns et al.l 1 19961: lEspanol & Reyengal l2003l : iLodato & Pried 


20 id) , with only I Watkins et al. ( 1996h suggesting that the two first 


derivatives approach is more accurate. By comparing our two im¬ 
plementations in Sec. ID we effectively compare both approaches. 
We find only small differences between the two approaches in terms 
of the overall accuracy, with the main advantages being that the di¬ 
rect second derivatives approach is both faster and easier to imple¬ 
ment. 


^ While IClearv & Monaghai] (Il999h proposed the harmonic mean, there 
is no detailed comparison between the two choices in their paper and the 
only proof that the harmonic mean correctly represents the second deriva¬ 
tive, apart from the numerical tests in their paper, involves a Taylor-series 
approximation where the harmonic mean reduces to the arithmetic mean. 


4 NUMERICAL TESTS 

A key issue in developing numerical codes for dust-gas mixtures is 
that there are few simple test problems that can be used to bench¬ 
mark the algorithm. We have partially resolved this issue by de- 
riving the analytic s olution for linear waves in such a mixture 
(iLaibe & Pried |201]]) and showing that the solution for a shock 
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in the limit where Av ^ 0 is the same as for the hydrody¬ 
namic case but with a modified sound speed (iLaibe & Pricell2012a| : 
iMiura & Glasslll982h . The dustybox solution lLaibe & Pricel201l[) 
is not relevant to this paper since we have already assumed that 
Av has reached its asymptotic value by using the terminal veloc¬ 
ity approximation. Hence, we use the dustywave and dustyshock 
problems to benchmark our algorithm. Our exploration of the dif¬ 
fusion approximation for dust suggested a new test problem with a 
simple analytic solution, which we describe in Sec. l4.3l 


4.1 Dustywave 

In DUSTYWAVE problem, we solve for the propagation of a linear 
wave in a dust-gas mixture . We set up the problem in ID as in 
our previous papers (lLP12al:lLP12bl:lLP14bh . using pd,o = Pg,o = 1 
(i.e. po - 2 and eq - 0.5) with a sinusoidal perturbation to the ve¬ 
locity and density of the mixture particles v(x) = vq sin(2:7rA) and 
p(x) = po [1 + ^Po sin(2;rA)], with amplitude vq = Spo = 1 X 10""^, 
with a corresponding thermal energy perturbation given by 6u = 
Po/Pgo^Pg- adiabatic equation of state is used with y = 5/3 and 
the thermal energy is set so that the initial sound speed Cs = 1. We 
use 100 SPH particles in the domain a E [0,1]. 

There is a fundamental inconsistency in the dustywave initial 
conditions when using the terminal velocity approximation because 
the setup of the problem and hence the analytic solution assumes 
that Avo = 0. By definition in the terminal velocity approximation 
we have Av = tsAf which is non-zero. Hence the solution even 
at t = 0 is not identical to the full one fiuid case. However, these 
differences become smaller at large drag and at later times. 

The numerical solution is shown after 4.5 wave periods in 
Figure [T] showing the g as and d ust velocities (filled and open cir¬ 
cles, respectively). As in lLP14bl we have reconstructed the gas and 
dust velocities on each particle from the barycentric variables, i.e. 
Vg = V - eAv and Vd = v -h (1 - 6)Av. The left Figure shows the 
results using the two first derivatives approach (Sec. l3.1l) while the 
right Figure shows the results using the direct second derivatives 
version (Se c. 13.41) in each case compared to the linear analytic so¬ 
lution from I lP 111 . There is no distinguishable difference between 
the two approaches. The solution in the regime where the terminal 
velocity approximation is valid (^ > 42; lower two panels in each 
Figure, corresponding to 4 > A^cour = 0.01) is within a few percent 
of the analytic solution. There is a conspicuous phase error at lower 
drag (^ = 10 and K = 1; first and second row), in part caused by 
the inconsistency in the initial conditions, which becomes worse as 
4 becomes larger, and in part because this is where the terminal 
velocity approximation breaks down. Nevertheless the general be¬ 
haviour in terms of the damping of the wave at intermediate drag 
is captured despite the inapplicability of the approximation in this 
regime. The behaviour at even lower drag (^ < 1; not shown) is 
incorrect; here the wave remains damped when using the terminal 
velocity approximation whereas the damping should decrease as 
the coupling tends to zero. Henc e the ful l one fiuid approximation 
should be used in this regime (e.g. iLPUhh . as we argued in Sec. l2.2l 

Figure [2| shows the solution for the density perturbation. Im¬ 
portantly, the analytic solution for density in the dustywave prob¬ 
lem quickly becomes nonlinear, particularly when the drag is weak. 
This can be seen by considering the limit of no drag: Assuming the 
dust is not submitted to any external force we have 

v(t) = v{t = 0) = Vo sin {kxo ), (46) 

implying 

x{t) - xq-vq sin {kxo) t. (47) 


Hence, from mass conservation, the dust density is given by 


Pd(0 ^ 


PdO (-^o) 


dxo 


PdO (-^o) 

|1 - vokcos (kxo) t\ ’ 


( 48 ) 


This result is physically consistent with the initial velocity profile: 
grains are depleted at x = ±7t, pile up at a = 0 and maintain a 
constant density at v = ±7t/2 (zero net flux of particles). In par¬ 
ticular, density fluctuations become of the order of the background 
on a typical ti me (vn ^)~^ and the analytic solution of the dustywave 
problem from I lP 111 cannot be applied anymore. It should be noted 
however that the velocities remain small and still agree with the 
solution of the linear problem. Hence, we have computed the refer¬ 
ence solution in Fig.|2]using a high-resolution (5 000 part icle) sim¬ 
ulation with our fully general one fluid al gorithm (lLP14bh : whereas 
in Fig. [T] we used the linear solution from lLPlH (both methods pro¬ 
duce indistinguishable results for the velocity field). 

Since there is no inconsistency in the density in the initial con¬ 
ditions, the solution using the diffusion approximation is more ac¬ 
curate for the densities than for the velocities (L 2 error of 0.06 at 
K - 100 and 0.006 at TT = 1000), though still becomes inaccurate 
(L 2 error > 0.5) for K < 10. As with the velocities, there is no dif¬ 
ference between the two implementations (compare left and right 
panels in Figure |2]), indicating that any inaccuracies are due to the 
physical approximation rather than the numerical scheme itself. 


4.2 Dustyshock 

The DUSTYSHOCK problem at strong drag was one of the most diffi¬ 
cult problems to solve using a two fluid approach due to the resolu¬ 
tion requireme nt h < that leads to overdampin g of the solution 
if not satisfied (lLP12ah . We have already shown in lLP14bl that this 
spatial resolution requirement is unnecessary when using a general 
one fluid formulation, although the drag still imposes a prohibitive 
timestep constraint, meaning that an implicit timestepping scheme 
(albeit a fairly simple one) is still necessary. Figure[3]show that with 
our present method we can capture the high-drag dustyshock so¬ 
lution using explicit timestepping without any timestep constraint 
other than the usual Courant condition, and without any particular 
spatial resolution requirements. 

W e set up the problem as usual, following the standard ISodI 
(Il978h shock tube with conditions in the gas for a < 0 given by 
(Pg,Vg,P) = (1.0,0.0,1.0) and for a > 0 given by (pg,Vg,P) = 
(0.125,0.0,0.125). We assume a constant dust fraction in the initial 
conditions (cq = 1), using 569 particles (corresponding to a particle 
spacing of Av = 0.001 for v < 0, an adiabatic equation of state 
with y = 5/3 and a drag coefficient K = 1000. In this regime 
the solution corresponds to the usual hydrodynamic solution with 
a modified sound speed (red lines in Figure [3]). In this respect we 
are testing only the ability of t he algor ithm to recover the zeroth 
order effect of a heavy fiuid fsee lLP14ar) . which from the results in 
Figure [ 3 ] can be seen to be true. 

As previously there is very little difference between the two 
implementations (comparing left and right panels) except that 
the direct second derivatives approach (c.f. Section [3Al ) produces 
slightly less noise in the pressure profile across the contact dis conti¬ 
nuity . While this can be important for some problems (see e.g. lPricel 
l2008h . such a minor difference is not enough to prefer this discreti¬ 
sation over the two first derivatives approach. However, given that 
the direct second derivatives algorithm is also significantly faster it 
may be preferred on this basis. 
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Figure 1. Gas and dust velocities (filled and open circles, respectively) in the dustywave problem using 100 SPH particles and our two implementations of 
the du st diffusion approximation: Two first derivatives (left) and direct second derivatives (right). These may be compared to the analytic linear solution from 
IlPIII given by the red solid (gas) and dashed (dust) lines. The L 2 error is within 6% of the analytic solution for ^ = 100 and within 2% for K = 1000, where 
the diffusion approximation is applicable (here for K > 42 corresponding to 4 > Afcour = 0.012). The solution becomes inaccurate at weaker drag (4 > 0.012). 
There is no discernible difference between the two implementations, except that the implementation with direct second derivatives (right) is faster. 



k _I_I_I_I_^_I_I_I_I_ J 

0 0.5 1 



k _I_I_I_I_^_I_I_I_I_ J 

0 0.5 1 


X 


X 


Figure 2. As in Fig.[T]but showing the density perturbation. The solution in this case ma y be com pared to the red solid (gas) and dashed (dust) lines, showing a 
high resolution non-linear solution computed using the general one fluid algorithm from lLPl^ . The solution is captured with increasing accuracy as the drag 
becomes stronger, with an L 2 error of 6% for ^ = 100 and 0.6% for K = 1000, but as expected becomes inaccurate in the regime where the approximation 
breaks down (4 > 0.012). 


4.3 Dustydiffusion 


Based on Equation [TSl we present a new test for dust-gas mixtures 
with a simple analytic solution. This consists of the steady diffusion 
of an overconcentration of dust. To set up the problem we consider 
a uniform density box with p = po = 1 and an isothermal equation 
of state P = CgPg with Cs = 1. In this case the dust diffusion can be 
described by Equation (TS] Eor the diffusion parameter we assume 
that the stopping time is a constant (this is equivalent to assuming 
an Epstein-like drag where 4 = p grainsgrain/ (pc s) is constant). 


4.3.1 Analytic solution 


The exact solution can be obtained by solving the equation 

^=V.(e; 7 Ve), (49) 

at 

where 77 = 4 Cg is a constant. We solve this by assuming spherical 
symmetry, i.e. 


de 77 d / 2 d 6 
dt dr \ dr 


(50) 


for which there are several known analytic solutions, including the 
general time-dependent solution 


6 (r ,0 =A\l 0 f]t + B\-^5 


r^ 

I0f]t + 5 ’ 


(51) 
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Figure 3. Results of the dustyshock test with a large drag coefficient, K=1000, comparing the use of two first derivatives to compute the dust diffusion (left) 
with the direct second derivative discretisation of the diffusion term (right). In both cases the numerical solutions agree with the analytic solution valid in 
the limit of infinite drag (solid red line), although the pressure is smoot her acros s the contact discontinuity when the diffusion term is computed directly. The 
advantage of the present scheme compared to the full one fluid ap proach |LP14bh is that we have used explicit timestepping. We also avoid the punitive h <t^c^ 
resolution requirement associated with the two fluid formulation (lLP12ah . 569 SPH particles were used. 



Figure 4. Dust fraction as a function of spherical radius in the 3D dust diffusion test at t = 0.0, 0.1, 0.3, 1, 3 and 10 (top to bottom) from simulations using 
50 X 58 X 60 particles. The numerical solution, projecting all particles in r, is given by the black dots and may be compared to the analytic solution given by 
the red lines. The left panel shows the solution with two first derivatives, while the right panel uses the direct second derivative. 


where A and B are arbitrary constants. We use this solution to verify 
our numerical scheme by solving only the diffusion equation via ei¬ 
ther Equations |^and[30]or Equation|43l with the particle positions 
fixed (for this problem onlyfl 


4.3.2 Results 



We set up the problem in 3D with 50 x 58 x 60 particles set on 
a uniform close-packed lattice in the domain v,y,z G [-0.5,0.5]. 
The positions of the y and z boundaries are adjusted slightly to 
ensure periodicity of the lattice across the boundary (the particle 

^ We attempted to construct an equilibrium situation involving all of Equa¬ 
tions [ToWT^ for example a hydrostatic equilibrium in a fixed potential. 
However, it is difficult to construct an equilibrium where the dust simply 
diffuses according to because the change to e causes a change to the 
pressure gradient and hence causes an acceleration to the barycentre also. 


Figure 5. Cross section of the dust density in the z = 0 plane in the 3D dust 
diffusion test at r = 0, 1 and 10 (left to right). 

spacing in v, y and z is Ap, V3A;?/2 and ^^6Ap|3 respectively, 
where Ap - 0.02). We use an isothermal equation of state, setting 
Cs = 1 and 4 = 0.1 such that 77 = 0 . 1 , and set the initial dust fraction 
using 
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Figure 6. Convergence in the dust diffusion problem, showing L 2 error for 
the solution within r < 0.2 as a function of the particle spacing. While 
both methods show second order convergence, the direct second-derivatives 
solution is more accurate because of oscillations in the two first derivatives 
approach propagating from the ‘kink’ in the initial € prohle seen in Figured 


consistent with Equation [5T] with B = 60 /r^ and A = €oBi. We set 
60 = 0.1 and Pc = 0.25. 

Figure m compares the numerical solution to the analytic solu¬ 
tion, while Figure|5]illustrates the general behaviour of the solution. 
The solution with Equation|43](right panel of Fig.|4]) is excellent (L 2 
error < 5 X10“^ for r < 0.2), apart from the physical deviation from 
the self-similar solution due to the transition to constant rather than 
negative e at the outer radius. The solution with using Equations |26] 
and [30] (left panel) is also good, but shows some low amplitude 
oscillations that develop from the propagation of the ‘kink’ in the 
initial epsilon profile. These oscillations are worse at lower resolu¬ 
tion (they can be smoothed out by adding some artificial dissipation 
in 6 but the solution is still not as good as using Equation l43l). 

Figure [ 6 ] quantifies these results with a convergence study us¬ 
ing 8 ^, 16^, 32^ or 64^ particles arran ged on a cu bic lattice. We 
show the L 2 error computed by splash ([Pricell200% from particles 
with r < 0.2. While in both cases the convergence is second order 
oc (bv)^, it can be seen that the direct second derivatives approach 
gives results more accurate by a factor of 5 at any given resolu¬ 
tion. Our results with both schemes when employing v instead of 
6 (Appendix |B]) are worse by a factor of ~ 2, again with a similar 
preference for the direct second derivatives approach. Thus while 
it is clear that all of our proposed numerical schemes correctly dis- 
cretise the diffusion equation, we find the discretisation using Equa¬ 
tion |43]to be more accurate for this problem. 


4.4 Dust settling in a protoplanetary disc 

Our final test is drawn from our intended application, namely the 
dynamics of small grains in protoplanetary discs. 


4.4.1 Setup 


We simplify the problem by considering only the vertical settling 
of grains in the r-z plane. That is, we set up particles in a two di¬ 
mensional cartesian box with an acceleration in the ‘vertical’ (z) 
direction given by 


GM 


(53) 


where we assume code units such that GM - 1 and set = 5 as 
a constant. The boundary conditions are periodic in the horizontal 
(x) direction and free in the vertical direction. We use an isothermal 
equation of state P - c^p where the sound speed Cs is set such that 

the aspect ratio H/Rq = /(Qq^o) = 0-05, where Qq = ^GMIR?^. 

The orbital time is therefore forb = 27t/Qo ^ 70 in code units. We 
set particles of equal mass initially on a uniform hexagonal lattice 
in the domain v G [-0.25,0.25] and z G [-3H,3H]. We specify 
the particle separation in the v direction to be either 16, 32, or 64, 
resulting in 16 X 56 = 856 particles at the lowest resolution, 32 x 
111 = 3552 particles at medium resolution and 64 x 222 = 14208 
at the highest resolution. 

We then stretch the particle distribution to matc h the equilib- 
rium density profile using the method described in [Price! (l2004l) 
where the z position of each particle is determined by solving the 
root finding problem 


f(z) = 


M(z) 

M{Zmax) 


(^0 ^min) 
^max “ -^min 


= 0 , 


(54) 


where M(z) = T p(z')dz', zo is the initial position of the particle 
and we set 


p(z)=poexp[-zV(2//2)]. 


(55) 


We set the mass of each particle equal to M(zmax) divided by the 
number of particles in the domain, consistent with the desired den¬ 
sity profi le. Equation |55] is a slight approximation (fourth order in 
z/7/; e.g. iLaibe et al.ll201^ but this is unimportant since we relax 
the particles into a hydrostatic equilibrium anyway, as described 
below. 

We set up the simulation initially with only gas and run the 
calculation to t = 1000 in code units (i.e. ~14 orbits) with both 
artificial viscosity and an artificial damping term of the form 





(56) 


where /damp = 0.03 in order to allow the distribution to relax to 
equilibrium. We then add dust to the simulation, assuming a dust- 
to-gas ratio of Pd/Pg = 0-01 by setting the dust fraction using 


Pd _ Pd/Pg 

p (1-fpd/Pg)' 


(57) 


We then evolve the simulation for a further 50-100 orbits. 

To give the problem physical meaning we consider a distance 
unit of lOAU (such that Rq = 5 corresponds to 50AU), a mass unit 
of IMq and the time unit set such that G = 1 in code units. This 
implies an orbital time of IttIQq = 353 years. A midplane density 
po of 10“^ in code units then corresponds to ^ 6 X 10"^^ g/cm^, giv¬ 
ing a disc surface density S ^ 55 g/cm^. We adopt a linear Epstein 
drag prescription, defining the stopping time according to Eq. [20] 
We set the intrinsic grain density Pgrain = 3 g/cm^. The midplane 
stopping time at Rq is given by 


tsQo = 1.35 X 10“ 


_3 / ‘pgrain \ / Pgrain 


; / “pgrain \ 

\ 1 mm/ 


3g/cm^ 


(58) 


4.4.2 Settling of millimetre grains 

We first perform a series of tests with a grain size of pgrain = 1 mm, 
chosen as a balance between the regime where the diffusion method 
is applicable and where it is still possible to obtain a solution in a 
reasonable time with the two fluid method. For the setup above 
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Figure 7. Settling of mm dust grains in a 2D (r-z) vertical section of a protoplanetary disc at Rq = 50AU (assuming H/R = 0.05; so Hq = 2.5AU) using 32x 111 
mixture particles. The plot shows dust density as a function of time. The top row shows the results using our new dust diffusion method. Th e solution may 
be compa red to that obtained with the full one fluid formulation from iLaibe & Pricel 12014bh (middle row) and with the two fluid formulation iLaibe & Pricel 
i2012dlM) (bottom row; uses 32 x 111 particles in both gas and dust). Our new method requires half the number of particles compared to the two fluid approach 
and is 50 times faster. 


this is at the limit of where the diffusion method is applicable, 
and indeed we found that Eq. [16] controlled the timestep, indicat¬ 
ing that the time depende nce of t he differential velocity has started 
to become important. In lLP12al we showed that it was necessary 
to satisfy h < to avoid overestimating the drag. For 1mm 
grains our two fluid calculations violate this criterion by a factor 
of ~ 9, 4.5 and 2.25 at the midplane at low, medium and high 


resolution, respectively, but do not appear to show overdamping. 
iLoren- Aguilar & Bate! (l2014l) found that the resolution problem is 
not as severe when the dust-to-gas ratio is low, suggesting that 
h < ecsts is a more precise resolution criterion. 

Figure [7] shows the dust density in the medium resolution cal¬ 
culations at intervals of ten orbital periods using three different 
methods. The top row shows the results with our new method em- 
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Figure 8. As in Figure |7]but showing the projection of dust density on all 
particles in 2D as a function of z, plotted at r = 0,10,20,30 and 40 orbits for 
the different methods. The direct second derivatives formulation (top) gives 
a slightly more smoothed result compared to the two first derivatives method 
(second row), the latter of which is indistiguishable from the solution with 
the full one fluid method (third row). The one fluid methods are less well 
resolved in the dust density for this problem but are give a significantly 
less noisy solution than obtained with the two fluid method (bottom). The 
diffusion method is 50 times faster at this resolution and grain size, with 
increasing performance gains for smaller grains. 


ploying the two first derivatives approach (Section 13.11) . The re¬ 
sults with direct second derivatives (Equation |43]) are similar but 
slightly less well resolved (see Fig. [U). The second row of Fig. |7] 
shows t he solution obtained with the general one fluid method from 
lLP14hl The main difference is that the differential velocity Av is 
explicitly evolved in that formulation, and so there is a timestep 
constraint from the stopping time which makes the simulation run 
~ 25 times slower for this grain size when computed with only ex¬ 
plicit timestepping with the constraint Ar < 4 . With the diffusion 
approximation the timestep constraint is inversely proportional to 
the stopping time although quadratically proportional to resolution 
(Equation [16]). The third row shows the solution obtained with the 
two fluid algorithm (lLP12ah . In this case instead of setting the dust 
fraction we added a separate set of dust particles copied from the 
gas particles but with 1% of the mass. This approach therefore re¬ 
quired twice the number of particles compared to our diffusion al¬ 
gorithm and the timestep is also constrained by the stopping time. 



Figure 9. Resolution study in ID version of dust settling problem (Figs.|7]& 
O, showing solution after 30 orbits using the diffusion approximation with 
direct second deri vatives ( top left), two first derivatives (top right), the one 
fluid method from lLPl^ (bottom left) and the two fluid method (bottom 
right). The particle number refers to the total number of particles in the 
domain in each case. In the one fluid methods the resolution follows the 
total mass rather than the dust mass, so the dust density is comparatively 
less well resolved. 


Computing this solution required approximately 50 times more cpu 
time than the diffusion method. 

All three methods produce dust settling on a comparable 
timescale, with only minor differences in the numerical solutions. 
A more detailed comparison is given in Figured showing the dust 
density on all particles as a function of z at the same times as those 
shown in Figure [T] The main noticeable difference is that the two 
fluid solution contains more noise in the particle distribution. This 
is because the dust is modelled as a separate set of particles that feel 
no mutual repulsion, compared to the one fluid case where the dust 
distribution benefits from the regular arrangement of the mixture 
particles. The approach with direct second derivatives (top row) 
produces a slightly over-smoothed solution compared to the the two 
first derivatives approach (second row) — with the latter giving re¬ 
sults that are indistinguishable from the full one fluid method (third 
row), showing that the diffusion approximation is indeed accurate 
in this regime. 

A major difference between the one fluid methods and the 
two fluid method is that resolution is tied to the total mass rather 
than the dust mass. This is evident in Figure [8] where the two fluid 
method (bottom row) can be seen to better capture the “wings” in 
the dust density at high latitudes. We quantify this further in Fig- 
urej^with a resolution study of the same problem performed in ID 
to avoid the particle noise in the two fluid approach. Settling means 
that after some time the dust covers a much smaller region of the 
domain than the gas, so the one fluid formulations under-resolve 
the dust compared to the two fluid method, since resolution is tied 
to the total mass, most of which remains at high latitudes. By con¬ 
trast in the two fluid approach resolution is tied to the dust mass and 
so is naturally placed towards regions of high dust density. This is 
both an advantage and a disadvantage to both types of approaches, 
it depends whether it is desirable to resolve the total mass or the 
dust mass. This point is discussed in lLP14bl mainly as an advan¬ 
tage to the one fluid method since it avoids the possibility of dust 
particles becoming ‘trapped’ below the resolution of the gas. 

In obtaining our results with the diffusion approximation, we 
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Figure 10. As in Figure|2]but comparing settling of grains of different sizes at in a vertical section of a protoplanetary disc at 50 AU. Each panel shows the dust 
density after 50 orbits computed with our new method (similar results are found with both implementations described in this paper). The decrease in settling 
time with increasing grain size (left to right) is clearly evident. Simulations in this regime are prohibitively slow with the two fluid approach. We used explicit 
timestepping in all cases, with the simulations for grain sizes of < 1mm constrained only by the Courant condition. 


found a few caveats to the numerical algorithm we derived in Sec¬ 
tion [3] First, we found it necessary with both to use the ker¬ 
nel for this problem to obtain a smooth and accurate solution (this 
kernel extends to 3h instead of 2h and so better approximates the 
Gaussian). To avoid particle pairing occurring at high latitudes with 
the quintic we set the ratio of smoothing length to particle spacing 
to 1.0 instead of 1.2, equivalent t o using a mean neighbour num¬ 
ber of 28.3 in 2D fsee |Pricell2012l) . The second caveat was that the 
dust fraction becomes negative around the edge of the collapsing 
dust layer in both of our implementations (and also with the full 
one fluid method). This arises because of the exact conservation 
of the dust mass in the algorithm, which causes a slight overshoot 
at the discontinuity. In order to smoothly handle this we derived 
an alternative approach (Appendix |B]) which guarantees a positive 
dust fraction, but we found it to give less accurate results than the 
method employing e (c.f. Sec. l4.3l) . Instead, we found that the most 
effective way of solving this was to simply set the dust fraction to 
zero on particles where it had become negative. This slightly vio¬ 
lates the exact conservation of the dust mass, but the error is small 
(~ 10"^ in 6 with the quintic) and it is a small price to pay for 
stability of the algorithm. 


4.4.3 Settling with different grain sizes 

Finally, we demonstrate the ability of the diffusion method to sim¬ 
ulate small grains in a protoplanetary disc by performing a series 
of calculations varying the grain size from O.l/zm to 1mm. Grain 
sizes below 1mm are difficult to simulate at all with the two fluid 
technique bec ause of the pun itive spatial and temporal resolution 
requirements ( LP12al : lLP12bh . With the general one fluid method 
presented in lLP14bl the limitation on the spatial resolution is re¬ 
moved, because we are no longer modelling the separation between 


fluids with physically separate resolution elements, but it is still 
necessary to use implicit timestepping. Yet these grains are impor¬ 
tant in protoplanetary discs as they control much of the thermal 
radiation. 

Figure [T0| shows the results of a series of medium resolution 
(32 X 111) calculations of dust settling for different grain sizes, 
shown after 50 orbits at Rq. We used only explicit timestepping, 
and for grain sizes smaller than 1mm the timestep was constrained 
only by the Courant condition. The different settling behaviour of 
the different grain populations in discs is clearly evident, with the 
micron and sub-micron grains remaining stuck to the gas at high 
latitudes, the millimetre grains settling effectively to the midplane 
and the 100 micron and 10 micron grains having partially settled. 

Being able to accurately and efficiently simulate single-size 
small grains in discs in this manner is the first step towards mod¬ 
elling an evolving grain population self-consistently. 


5 DISCUSSION AND CONCLUSIONS 

We have derived and implemented a numerical scheme for de¬ 
scribing the dynamics of small dust grains coupled to a gas, us¬ 
ing Smoothed Particle Hydrodynamics, in the limit where the stop¬ 
ping time is short compared to the computational timestep. This re¬ 
quires solving one additional diffusion equation as well as the usual 
equations of hydrodynamics slightly modified by some additional 
terms. We derived two implementations, one where the diffusion 
equation is computed using two first derivatives (Section [341 and 
one where direct second derivatives were employed (Section [34l . 
We found only minor differences between the two approaches on 
the test problems we tried. Given this, we recommend the direct 
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second derivatives approach (Sec. I3.4t . which is both simpler and 
faster because it does not require an extra loop over the particles. 

As discussed in Sec. 12.21 the terminal velocity approximation 
or, as we prefer, the “diffusion approximation for dust”, is valid 
when the stopping time is less than the computational timestep. The 
simple way to guarantee this validity in practice is to ensure that the 
diffusion timestep (Eq. [Tb]) is not constraining the timestep , other- 
wise the more general one fluid approach implemented in iLPldhl 
where the time dependence of Av is kept should be used instead. 
In this sense the method we h ave desc ribed is complementary to 
both the full one fluid approach iLPldbll and the two fluid approach 
(lLP12ah . The main difference is that the other methods need im¬ 
plicit timesteps when the grain size is small (4 < Ar), whereas 
this method requires implicit timesteps when the grain size is large 
(4 > AO, but this is where the approximation breaks down anyway. 

Finally, we considered only one grain size at a time in this pa¬ 
per. We have recently generalised our one fluid formulation to de¬ 
scribe an arbi trary nu mber of grain populations all within a single 
fluid mixture (lLP14ch . Our next step will be a numerical implemen¬ 
tation of this more general formulation, including the simplification 
to a diffusion approximation, as well as modelling the evolution of 
the grain population including growth and fragmentation. 
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APPENDIX A: PROOF THAT EQUATION|4l]IS A DISCRETE FORM OF EQUATIONH 


Here, we prove that the expression obtained for the second term in Eq.|45]by enforcing the conservation of energy, namely 


^ ou \r^u\ 


2(\ - e,)pa Pi “ "'Kiy 

is indeed a discrete form of the corresponding term in Eq.|4l i.e. 


-VP-Vu. 

Pg 

We proceed, following IPricd (l2012h . by identifying -2Fabl\rab\ as equivalent to the second derivative of a (new) kernel function, i.e. 


(Al) 


(A2) 


V^Yab ^ (A3) 

\rab\ 

It may be sho wn straightf orwardly that this new kernel Yab indeed satisfies the normalisation conditions appropriate to the kernel second 
derivative rsee lPricell2012l for more details). We can then take the Laplacian of the standard SPH summation interpolant with this kernel, i.e. 

(A4) 

V Pb 


to give 

b P” 

By writing dAlTi in the form 

A y - Ub){Da + Db)(P, - PbW^Y,b, 

4p“ ^ Pb 

we can then use (IA51 ) to translate the various terms. Expanding (IA6I) we have 

"j ~ ^ j {Pa^a^a ~ Pa^b^^a Pa^aP^b ~ Pa^bP^b ~ Pb^aP^a Pb^bP^a ~ Pb^aP^b Pb^bP^bli^ Pab- 

4Pg V 

Translating each of the terms in turn using ( lA5t gives 


A ^PuDV^l - PDV^u + PuV^D - PV^iuD) - uDV^P + DV^(Pu) - uV\PD) + V^(FmD)] . 


Expanding the V^(ab) terms using the vector identity 
V^(ab) = aV^b + 2(Va • Vb) + bV^a, 
and expanding the last term using 


V^(PuD) = uDl^P + PDl^u + Pul^D + 2w(VP • ID) + 2P)(VP • Vw) + 2P(VP) • Vw), 
we find, upon simplification that (IA81) reduces to simply 

A[4D(VF.Vm)]. 

4pg 

Hence, ( lA6t and so dSD is a discrete form of 

D eL 

— (VP -Vw)= —VP Vu. 

Pg Pg 

QED. 


(A5) 

(A6) 

(A7) 

(A8) 

(A9) 

(AlO) 

(All) 

(A12) 


APPENDIX B: ENFORCING POSITIVITY OF THE DUST FRACTION 

While the usual SPH density summation enforces positivity of the total density, in the one fluid approach there is no constraint on the 
positivity of the dust fraction, being simply evolved via a differential equation. We found during our testing of the algorithm (Section |4^ 
that this can occur in practice, even though we conserve the total dust mass. An simple example is where e is non-zero on only a fraction 
of the particles and zero on others, implying an infinite gradient in e at the discontinuity surface, which as the dust front evolves can lead 
to negative e on the particles that initially had zero. It should be noted that however that those errors are small and kernel dependant (i.e. of 
order 10"^ with a quintic kernel). We discuss other possible solutions in Section l4!4l but here present one such solution, which is to evolve 
the quantity 

s=^/^, (Bl) 

instead of 6. We can enforce the same conservation of dust mass but with a guaranteed positivity of the dust fraction since 

ea = sl/Pa. (B2) 
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B1 Continuum equation 

In terms of 5, the local equation for dust mass conservation is 


dt 


-£jv- 


f) 


Av 


^ 1 


-lAv 


Vv, 


-]av-V. 


(B3) 


where as previously d/df is the convective derivative using the barycentric velocity. For the case of hydrodynamics rSection [2.1.2l) where 
Av = 4VP/[(1 - lp)p], we can simplify this to 


d^ 

d7 


„ , st^VP\ t. „ 

V • I- H—VP • 

P j P 


- 5 V.V. 


(B4) 


B2 SPH implementation with two first derivatives 

Conservation of the total dust mass implies 


dt 

\ a 

giving 


^ —s 


Pa 

fUa d^, 


= 0 , 


\ 1 fUa 05 '^ \ ^ ^a ^Pa 2 


^ Pa dt 

To enforce Eq. IB6I we compute the evolution of according to 
difl Pa V / 1 “ ^a/Pa , „ ^ ^ “ ^/Pb 


Q^a Pa V 


iiaPl 


AVa • VaWab iK) + i-f^AVfc • V,Wa, (h,)] + V m, {Va - Vfc) • VaWab (ha) • 

^hpl I ^Pa^a V 


(B5) 


(B6) 


(B7) 


^ \ -ai-a ^bPij / -i-a—a ^ 

The first term of the right-hand side of Eg. IB7I corresponds to the first two terms (i.e. inside the brackets) of the right-hand side of Eg. lBSl in 
the continuous limit (note the factor inside the SPH summation). The contribution of this term to the left hand-side of Eg. IBbl is zero as it 
leads to a double summation of an antisymmetric term with respect to the indices a and b. The second term of the right-hand side of Eq. IB7I 
corresponds to the V • v term of Eg. lB3l and provides the right hand-side of Eg. IBbI which can be seen by differentiating the SPH density 
summation (US with respect to time. 


B3 SPH implementation with direct second derivatives 

We can also construct a method evolving 5 but with direct second derivatives. Here we discretise Eq. lB4l using 


dsg 

dt 


^- iPa + I^b){Pa ■ 

/, Pab 


p \ Sg 

^\rab\ 2pa^a 


'y ^b^ab ' ^^ab(ba)-> 


(B8) 


where D = s^ts/p as previously. It is straightforward to show that this is indeed a discretisation of (IB41) using the method described in 
Appendix I aI The average density is required in the denominator in order to conserve the total dust mass, which can be verified by substituting 
dM]) into ( lB6t . 
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